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Ground-state properties of trapped Bose-Einstein condensates: 
Extension of the Thomas-Fermi approximation 
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We derive general approximate formulas that provide with remarkable accuracy the ground-state prop- 
erties of any mean-field scalar Bose-Einstein condensate with short-range repulsive interatomic interac- 
tions, confined in arbitrary cylindrically symmetric harmonic traps. Our formulation is even applicable for 
condensates containing a multiply quantized axisymmetric vortex. We have checked the validity of our 
formulas by numerically solving the 3D Gross-Pitaevskii equation. 
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I. INTRODUCTION 

Since the experimental realization of the first dilute 
Bose-Einstein condensates (BECs) of trapped atomic gases 
JH 0, 111] there has been great interest in the study of the 
physical properties of ultracold quantum gases [4]. Under 
the usual experimental conditions, these systems can be ac- 
curately described by the Gross-Pitaevskii equation (GPE) 
JUl, a mean-field equation of motion governing the behav- 
ior of the condensate wave function ip(r). In the stationary 
case the GPE is given by the nonlinear Schrodinger equa- 
tion 



2m 



V 2 + V(r) + gN \ 



ip = flip, (1) 



where \x is the chemical potential, N is the number of 
atoms, g = An h 2 'a/ 'm is the interaction strength, deter- 
mined by the s-wave scattering length a, and V(r) is the 
trap potential. 

In this work we shall concentrate in the usual case of 
condensates with repulsive interatomic interactions (a > 
0), confined in cylindrically symmetric harmonic traps, 
V(r) = \m(u\r\ + uj 2 z z 2 ). 

To obtain the condensate ground-state properties one has 
to solve the nonlinear differential equation ([]]). No ex- 
plicit analytical solutions are known, so that, in general, 
this has to be done numerically. However, solving numeri- 
cally the 3D GPE is a nontrivial computational task, espe- 
cially for highly asymmetric trap geometries, where large 
basis or gridpoint sets can be required to guarantee conver- 
gence. Fortunately, approximate analytical solutions can 
be found in certain limiting cases. In the Thomas-Fermi 
(TF) regime, which essentially occurs for condensates with 
a large number of atoms, the kinetic energy can be ne- 
glected to a good approximation. In this case the GPE re- 
duces to a simple algebraic equation and one can obtain 
explicit analytical expressions for the condensate ground- 
state properties [6]. This approximation, however, cannot 
reproduce correctly the decay of the wave function at the 
boundary of the atomic cloud, a region where the kinetic 
energy has a decisive influence. Corrections to the TF ap- 
proximation have been proposed to account for the proper 



behavior of the density cloud at the condensate surface 
flULIsl]. In the (ideal gas) perturbative regime, when the 
number of atoms in the condensate is sufficiently small, ex- 
plicit analytical expressions can also be obtained by treat- 
ing the interaction energy term in Eq. (Q]) as a weak pertur- 
bation. 

Another approximation scheme that has proved to be 
useful in the characterization of dilute BECs is that based 
on variational techniques [6]. This approximation method 
has the additional interest that it can be e qually applied to 
the study of the condensate dynamics floL UlTl . By using 
appropriate variational trial wave functions, the ground- 
state properties of trapped Bose-condensed gases have 
been studied beyond the two analytically solvable regimes, 
for both isotropic lfl2ll and highly anisotropic condensates 
1 130. For this purpose semiclassical approximations (h — ► 
0) have also been considered lfl4l[T5ll . These approximate 
methods can be applied to the study of vortex states as well 
EQSQ1. 

Particular attention has been devoted to the physical 
prop erties of quasi- ID lfl9l 2fjl 2ll I22I l23ll and quasi-2D 



prop 

\ 24l 12511 condensates, systems so tightly confined in the 



radial or the axial dimension, respectively, that the corre- 
sponding dynamics is restricted to zero-point oscillations. 

In a previous work, by modifying the usual TF approx- 
imation conveniently, we derived very accurate approxi- 
mate analytical expressions for the ground-state properties 
of trapped spherical, cigar-shaped, and disk-shaped con- 
densates with an arbitrary number of atoms in the mean- 
field regime lEfjll . In this work we extend our previous 
results and derive general approximate formulas that pro- 
vide with remarkable accuracy the ground-state properties 
of any mean-field scalar Bose-Einstein condensate with 
short-range repulsive interatomic interactions, confined in 
arbitrary cylindrically symmetric harmonic traps, and even 
containing a multiply quantized axisymmetric vortex. 



H. MODEL 

We consider a BEC with an axisymmetric vortex line of 
topological charge q, and confined in a harmonic trap char- 



2 



acterized by oscillator lengths a± = ^h/mLO± and a z = 
\J h/muj z . The trap aspect ratio is given by A = u) z /ujj_. 
In this work we shall distinguish between ground and vor- 
tex states only through the value of the vortex charge. 
Accordingly, we shall refer to the lowest-energy state of 
the condensate compatible with an axisymmetric vortex of 
charge q as the ground-state with a charge-q vortex. This 
terminology will permit us to make a unified treatment in 
which the actual ground state of the condensate is simply a 
particular case corresponding to q = 0. 

For condensates with a sufficiently small number of 
atoms the mean-field interaction energy can be treated as a 
weak perturbation. In this perturbative regime the conden- 
sate wave function that minimizes the energy functional is 
given, to the lowest order, by 

ip q (r±,z, 9) = exp(iq0)tp q (r ± )(j>(z), (2) 

with 

<p 9 (r ± ) = (7ral\q\\)-Hr ± /a ± ) M eM-ri/2al), (3) 



0(z) = (^)- 1 / 4 e X p(-z72aj). (4) 

Using this result one obtains the following analytical ex- 
pression for the chemical potential: 

1 



-hw z + (|g| + l)hw± + gn, 



(5) 



where n = c q N/ (27r) 3 / 2 a\a z is the mean atom density 
and the coefficient c q is a function of the vortex charge that 
takes values in the interval [1, 0) and accounts for the di- 
lution effect that the centrifugal force associated with the 
vortex has on the condensate mean density 

q 2^\{\q\\y K) 

For q = one has c q = 1, and the usual results for the 
ground-state properties of a BEC in the perturbative regime 
are recovered. For a unit charge vortex the dilution coeffi- 
cient takes the value 1/2, and, in general, for |g| 3> 1, it 
decreases slowly as l/^/vrlgj. 

In the Thomas-Fermi regime, for condensates in the 
ground state (q = 0) and with a sufficiently large number 
of atoms, one can neglect the kinetic energy in comparison 
with the interaction energy, and the stationary GPE leads 
to 



1 1 

/i = -moj 2 z z 2 + -miv 2 ± r 2 ± + gN \i/j(r±, z)\ 



(7) 



In the presence of a vortex of charge q ^ the above 
equation is no longer a good approximation. In this case 
neglecting the kinetic energy amounts to neglecting the 
vortex itself. However, for a large N, the condensate den- 
sity cloud outside the vortex core is still given, to a very 
good approximation, by the Thomas-Fermi profile. We 
thus shall assume the above TF expression to be valid up to 
a lower cutoff radius 



r\ = V2(M + 1) 



determined from the condition that, in the presence of a 
vortex, the contribution from the radial harmonic oscillator 
energy cannot be smaller than (\q\ + l)fkv±, i.e., 



1 2 / \2 



(\q\+l)tiu ± . 



(9) 



Likewise, the contribution from the axial harmonic oscil- 
lator energy should not be smaller than the corresponding 
zero-point energy. To account for this fact we introduce a 
second cutoff z , defined through the condition 



1 



,2 „2 



1 



-hco. 



(10) 



which yields the axial cutoff 

z Q = a z . (11) 

This leads us to consider the TF expression (0 applicable 
in a volume V 3 that corresponds to the usual TF ellipsoidal 
density cloud truncated at both rj_ = r° and \z\ = Zq, 



+ z 2 /Z 2 F < 1, 



r i > r A \z\ 



(12) 



V 3 = {(r ± ,z): r 2 jR 2 

> z a } 

y / 2fi/hujj_ a x and 
are the usual TF values for the con- 



In the above equation i? TF 



densate radius and axial half-length, respectively. In our 
approach, however, the condensate radius and axial half- 
length do not coincide with the above TF expressions. Be- 
cause of the cutoffs we have introduced they are given in- 
stead by 



R = y / 2[i/hu A 



A a 



(13) 



(8) 



Z = ^2fi/hw z -2{\q\ + 1)/A a z . (14) 

However, when fi ^> (\q\ + l)huj±, then Z becomes indis- 
tinguishable from Z TF . This is true, in particular, when- 
ever A 3> 2(\q\ + 1). Likewise, if /i 3> \hio z , then 
R — > Rtf- 

More generally, for condensates with /x 3> + 
which occurs whenever A ^> 2(\q\ + 1), the relative con- 
tribution to the condensate properties coming from the re- 
gion r_L < r° becomes negligible. The same occurs with 
the relative contribution from the region corresponding to 
\z\ < z when the chemical potential is much larger than 
the axial zero-point energy, which occurs whenever A <C 1. 
It is then clear that when the number of atoms is sufficiently 
large that [i S> (\q\ + + \^>z the dominant contri- 

bution comes from the V3 volume introduced above. 

Only for condensates with fi ~ (\q\ + l)huj± is the 
contribution from the region r±_ < r° the most signifi- 
cant one. These are condensates with the transverse dy- 
namics frozen in the lowest energy state compatible with 
a charge-q axisymmetric vortex. In this case the conden- 
sate wave function can be factorized as ip q (rj_,z,8) = 
exp(iq6)if q (r ± )(p(z), with f q (r±) given by Eq. (O. Af- 
ter substituting this wave function in the stationary GPE, 
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multiplying by tp*, and integrating over the radial dynam- 
ics, one obtains 

^mu 2 z 2 + {\q\ + l)hw ± + g m N\<P{z)\ 2 = M) (15) 

where g 1D = c q g /2Tra 2 L , and we have neglected the axial 
kinetic energy (~ ^fau) z ) against (\q\ + l)hui±. It is conve- 
nient to rewrite g 1D as gK{ x n 2 , where n 2 = l/vr(r° ) 2 = 
l/[27raj_(|g| + 1)] is auniformmean density per unit area 
normalized to unity in the region r_|_ < r° , and 

k- 1 = (\q\ + l)c q (16) 

is the appropriate renormalization factor. The important 
point is that the effect of the vortex can be incorporated 
in an exact manner into a localized uniform mean density 
n 2 with a renormalized interaction strength gn^ 1 . We shall 
assume that, to a good approximation, this still remains true 
for condensates with an arbitrary chemical potential. 

On the other hand, the contribution from the region 
\z\ < Zq is the most significant one only for conden- 
sates with [i ~ In this case the axial dynam- 
ics is restricted to zero-point oscillations, and the wave 
function of such quasi-2D condensates can be written as 
ifj q (rj_,z,6) = exp(iq6)cp(r±)(f)(z), where now 4>{z) is 
given by Eq. ©. Substituting into the stationary GPE and 
integrating out the axial dynamics one obtains 



1 



'-fwj. 



1 



mulr 2 x + g 2D N\(p(r ± )\< 



Mi 



(17) 



where g 2 n = g/v2na z , and now we have neglected the 
radial kinetic energy against the axial zero-point energy, 
which is a good approximation as long as \i ~ \hw z . As 
before, it is convenient to rewrite <? 2 d m terms of a uni- 
form mean density per unit length normalized to unity in 
the volume \z\ < z . To this end we introduce a renormal- 
ization factor k 2 x = y2/7r and rewrite g 2 v as gn 2 x n\ 
with fix = l/2a z . This indicates that the contribution 
from the axial zero-point oscillations, which is the dom- 
inant contribution in these quasi-2D condensates, can be 
properly accounted for by simply introducing a localized 
uniform mean density per unit length with a renormalized 
interaction strength. Again, we shall assume this to be also 
valid for condensates with an arbitrary chemical potential. 
One finds, however, that somewhat more accurate results 
are obtained when one lets the renormalization factor k 2 1 



to approach unity in the TF regime 12611 . Since the final 
results are little sensitive to the specific functional form of 
k 2 x , we propose one of the simplest possibilities JIT 



Rtf{x? 



, (18) 



where Q(x) is the step function, \2 = Na/X 2 a z , and 
Rtf(X2) = (15x2) 1 ^ 5ffl ± is the TF radius. 



Motivated by the above ideas and the fact that there ex- 
ists a direct relation between the number of atoms and the 
size of a trapped BEC, we propose the following ansatz for 
the ground-state properties of any mean-field scalar Bose- 
Einstein condensate with an axisymmetric vortex of charge 
q, confined in an arbitrary axisymmetric harmonic trap: 

r G V 3 
r£V 2 
r G Vi 
reV 

with ip = elsewhere. In the above equations a q = (\q\ + 
1) and n = K1K2, with Ki and k 2 defined by Eqs. ( fT6l l 
and ( fT8l . respectively. As already seen, fh\ = l/2a z and 
n 2 = l/2Tra 2 L a q , whereas n is an effective mean density 
(per unit volume) localized in V and defined in terms of 
\i through the above expressions. Note that Kq 1 is exactly 
the renormalization constant required to make the latter of 
the above equations compatible with the other ones as well 
as with the perturbative result (0. The outer volume V 3 is 
defined by Eq. (fT2l while the remaining inner volumes are 
defined by 

V 2 = {(r ± ,z): r° < r± < R A \z\ < z } , 



\muj 2 z z 2 + 
2 


^mu} 2 ± r 2 ± + gN \ip(r ± , z)\ 2 


= V 


-huj z + ^muj 2 ± r 2 ± + gn 2 1 Nn 1 \ip(r x )\ 2 


= V 


-muj 2 z 2 + 
2 


a q hv±_ + gn x 1 Nn 2 \(j)(z)\ 2 


= V 




^hw z + a q hw ± + gK~ l nv 


= V 



z < 



V x = {(r x ,z): r ± <r° L A z < 

V = {(r x ,z): r^<r° ± A \z\ < z } . 

The ansatz we have just introduced represents a direct 
generalization of our previous proposal in Ref. Il26ll and 
extends the applicability of the approach to mean-field con- 
densates confined in axisymmetric harmonic traps with an 
arbitrary geometry, and containing an axisymmetric vortex 
of charge q. As we shall see, in the appropriate limits the 
results obtained in the present work reduce to those previ- 
ously obtained in Ref. [26]. 

From Eq. (fT4l l the chemical potential can be written in 
terms of the dimensionless axial half-length Z = Z/a z as 

fj, 1—2 1 



r,, \* + a (I ' 



+ 1 



(19) 



As usual, the condition that the condensate contains N par- 
ticles determines the precise value of /x. After a straight- 
forward calculation one obtains 

15 

1 / _ - 1 zr 

2 \ A 2 2 



Xo 

\5/3 



z + 

Ml 



6 



,2 1 



A 



4 



(20) 



where (3 q = Kx(\q\ + 1) = l/c q , £ n = (« 2 - 1/n) with 
n = 1, 3, 5, . . ., and we have defined the dimensionless 
interaction parameter Xo as 

Xa = Na/a , (21) 



4 



with a = (a^a z ) 1 / 3 being the mean oscillator length. 
From Eqs. ( PT3l) and (fT4l one also finds the following ex- 
pression for the condensate radius R = R/a±: 



R 



X(Z -l) + 2(|g| + l). 



(22) 



The mean-field interaction energy per particle e int = 
E int /N is defined by e int = (1/2) / g r N |^(r)| 4 d 3 r, 

where N \ip\ 2 represents the local density in each region 
and g r denotes the corresponding renormalized interaction 
strength. After some calculation one obtains 

A 5 /3 f 8 -7 £ l7 =6 



£int 



8Xo 



105 



-Z 



2±Z 
6 



+ 



8/3,-^5 



+ 





£ 3 


A 


2 




£7 


A 


6 



^ z 



15A 
A 



k 

4 



- Z 



(23) 



Finally, the kinetic and potential energies can be readily 
obtained in terms of the previous results by using the exact 
relations |0] 

e ki n = E kin /N = fi/2 - (7/A)E int /N, (24a) 
e P ot = E pot /N = n/2 - (l/4)E- mt /N. (24b) 

Equations ([T9ll- (l24]| provide the ground-state properties 
we are looking for. All that is needed is to solve the quintic 
polynomial equation (l20l l. This is a general equation that 
provides the axial half-length Z of any mean-field scalar 
condensate as a function of only three parameters: the in- 
teraction parameter Xo> the trap aspect ratio A, and the vor- 
tex charge q. In certain particular cases it is possible to find 
useful approximate analytical solutions. However, in gen- 
eral, Eq. d20T ) has to be solved numerically. It is important 
to note that this is a trivial computational task that can be 
done immediately with the built-in capabilities of symbolic 
computational software packages such as MATHEMATICA 
or MATLAB. In fact, to obtain the roots of a polynomial one 
simply has to type in a single instruction and the answer is 
instantaneous. 



III. LIMITING CASES 

The above formulas simplify considerably in two limit- 
ing cases that, essentially, correspond to condensates con- 
fined in disk-shaped traps satisfying A ^> 2(\q\ + 1) 
and cigar-shaped traps satisfying A -C 1. We have al- 
ready found these two limiting cases before. As men- 
tioned above, in the first case the relative contribution to 
the condensate properties coming from the inner cylinder 
r± < r° becomes negligible, while, in the second case, it 
is the relative contribution from the inner disk \z\ < z that 
becomes negligible. Under these circumstances we shall be 
able to find approximate analytical solutions of the polyno- 
mial equation (l20l l. 



Disk-shaped traps 

Taking into account that j3 q /\ < {\q\ + 1)/ A, it follows 
that in the limit A > 2(\q\ + 1) Eq. reduces to 



y 2 = — Z H Z Z H . 

A 15 8 4 8 



(25) 



where X2 = Na/\ 2 a z is now the only relevant physical 
parameter. Using Eq. (l22l one can easily see that for q = 
the above equation coincides exactly with that obtained 
previously in Ref. Il26ll . This is true in general (i.e., for any 
q) whenever A 3> 2(\q\ + 1). Under these circumstances 
the contribution of the vortex can be neglected to a good 
approximation and we can use the analytical solution found 
in Ref. JH 



Z 



1 + 



(1/15 X2 ) 8/5 + (k 2 /8 X 2Y 



-1/4 



(26) 



From this result one immediately obtains the chemical po- 
tential using Eq. (Q3D> 



JL 



¥■ 



(27) 



On the other hand, in the limit we are considering, the 
interaction energy d23l) becomes 



€int 



z 7 + ^-z 6 

8x2 v 105 " 6 



+ kr 
2 2 



6 



(28) 

As before, it is not hard to see that this equation is the same 
as that obtained previously in Ref. J2(J. 

Another relevant physical quantity in the characteriza- 
tion of disk-shaped condensates is the condensate density 
per unit area, defined as n 2 (rjj = N J dz \t/j(r±, z)\ 2 . 
Outside the vortex core (r± > r"), which we are neglect- 
ing in this limit, a straightforward calculation leads to 

„ i(r±)= 6feK)-i] + [2M,K)r-i i (29) 



4naa, 



6rraa, 



where n 2 (r± > R) = and fi z {r±) = /j. z (r ± )/hio z is 
given by 



^)^ + l(R/Vx) 2 (l-^ r 



(30) 



— 2 
with R ~ 



\(z 2 



1) [Eq. d22l l. Again, this result co- 
incides with that derived in Ref. [26]. In fact, in general, 
in the limit A S> 2(\q\ + 1) the formalism proposed in the 
previous section becomes independent of both A and q and 
reduces to that developed in Ref. [26]. 

It can be easily verify that fi z (rj_) is nothing but the lo- 
cal chemical potential, defined as fj, z (r± ) = fi— 

lfl5ll . Equation (l29l l is a cubic equation in Jl 1 / 2 which has 
only one real solution. Solving this equation, after some 
algebra one finds the following expression for the local 



5 



chemical potential as a function of the condensate density 
per unit area ll27ll : 



1 



v 



(31) 

where ij = 4+6£i — £f + 24:Traa z n 2 (r±). In the TF regime 
(X2 3> 1 ^ aa z n 2 3> 1), the above equation reduces to 



2/3 



(32) 



which coincides with the expression that can be obtained 
directly from the 3D Gross-Pitaevskii equation in this 
regime. In the quasi-2D perturbative limit (% 2 "C 1 — > 
aa z n 2 <C 1), Eq. (l3ll reduces to 



(33) 



H z (r±) = 1/2 + 2v27raa z n 2 (r 



This is again the correct result, as follows from the pertur- 
bative solution of the GPE in this limit. 

Equations (|29ll-(l3ll permit us to derive a general for- 
mula for the (local) radial (first) sound velocity c 2 d of a 
disk-shaped condensate, which is defined by 



n 2 dfi z 



'2D 



From Eq. 



mc 



2D 



m dn 2 
one immediately obtains 

6 [2&(r x ) - 1] + [2/Z,(r x )] 



3/2 



36+3[27Z 2 (r ± )] 



1/2 



(34) 



-■ (35) 



This equation in the TF regime (% 2 S> 1) reduces to 



mc 



2D 



2 -r ^ 



2^ 

-^=aa z n 2 (r ± ) 



2/3 



(36) 



while in the quasi-2D perturbative regime (% 2 "C 1), it 
reduces to 



"T-C 2D _ 1 

-^7 = ^ ) -2 



2V27raa z n 2 (r ± ). (37) 



These are the correct limits as follows from the substitution 
of Eqs. (1321 and (l33l l into Eq. d34l . In fact, it is not hard to 
see that all the analytical formulas derived in this section 
reduce to the correct expressions in both the TF and the 
perturbative regimes 12611 . 

For the particular case of a homogeneous disk-shaped 
condensate (no radial confinement and n 2 constant) the 
authors of Ref. 12311 obtained an expression for the ra- 
dial sound velocity that interpolates between the above two 
regimes. Such expression leads to the correct perturbative 
result and reproduces to a good approximation the TF re- 
sult. Even though this expression is analytic, it is too com- 
plicated to be written in a useful compact way and the au- 
thors of Ref. II23I1 do not provide an explicit analytical for- 
mula in their work. 



Cigar-shaped traps 



In the A < 1 limit Eq 

Xi 



reduces to 



±-(j\zT + l -p q {j\zy 



(38) 



with Xi — XNa/a±. In this limit the formalism becomes 
independent of A and the vortex contribution enters in a 
rather simple way through the parameter 



Pq 



2 2 l«l(|g|!) 
(2M)! 



(39) 



This fact will permit us to find an approximate analytical 
solution. 

Given a polynomial equation P(x) = x> we define the 
residual error associated with the approximate solution x £ 
as (P(x £ ) — x)/x L26]. We have explicitly verified that 
the expression 



1 



+ 



1 



+ 



1 



(15 X i) f + § 57X1 + 345 (3 Xl /& 



(40) 

satisfies Eq. d38l) with a residual error smaller than 3.2% 
for any Yi G [0, 00) and < \q\ < 10. In fact, as seen 
in Ref. 112611 . in the absence of vortices (q = 0) the above 
solution is somewhat more accurate and the error is less 
than 0.75%. 

From Eq. ( fT9l l the chemical potential is given now by 



{\ q \ + l) + l -(^\Zf 



In the A <C 1 limit the interaction energy becomes 



£int 



1 



15X1 



U^fxzy + p q (Vxzf 



(41) 



(42) 



On the other hand, the condensate density per unit length, 

(r±,zy 2 * 



ni(z) = N J 2irrj_drj_ 



4a 



z' 
~Z? 



+ 



, is 

(VAZ) 4 
16a 



with rii(z) = for \z\ > Z. 

The local chemical potential, defined as H±_{z) 
\muj 2 z z 2 UM, is given by 



z" 

~Z? 
(43) 



= fi 



(\ q \ + i) + ±(V\zy 



1 



z 

'z 2 



(44) 



The above analytical formulas generalize those obtained in 
Ref. l26h to the case of condensates containing an axisym- 
metric vortex. This is particularly interesting because, in 
this case, the usual TF approximation does not lead to ex- 
plicit analytical formulas for the condensate properties. In 
the absence of vortices (3 q — > 1 and one recovers the results 
of Ref. IH. 



6 



Substituting Eq. (l44l into Eq. d43l) and solving for 
~p>x.{z) = fj>±(z)/hujj_ one obtains the local chemical po- 
tential as a function of the condensate density per unit 
length 



£i» = (M + 1) + \//3 9 2 + 4am (» - A 



(45) 



This equation in the absence of vortices takes the simple 
form 



In the TF regime (xi S> 1 
sion reduces to 



x /l + 4an 1 (z). (46) 
an! 3> 1) the above expres- 



Hj_(z) = 2^/anJz), 



(47) 



which is the well-known result that can obtained directly 
from the GPE in this regime [22]. In the mean-field quasi- 
1D limit (xi <C 1 -> arai < 1) Eq. (06]) reads 

Ji ± {z) = l + 2an 1 {z). (48) 

This is again the correct result that follows from the pertur- 
bative solution of the GPE in this limit 12211 . In fact, it can 
be easily verify that all the analytical formulas derived in 
this section have the correct limits in both the TF and the 
perturbative regimes l26Tl . 

From the above results one can derive a general formula 
for the (local) axial (first) sound velocity Cid of a cigar- 
shaped condensate, defined by 

2 ™i <9m_l 



'ID 



m drii 

Substitution of Eq. (03]) in Eq. (09]) leads to 



(49) 



mc 



ID 



I 4a 2 nl(z) 
/3. 2 + 4an 1 (z)' 



(50) 



In the absence of vortices and in the TF regime (/3 q 
Xi 1) the expression above reduces to 



1, 



mc 



ID 



1 



(51) 



This result is in agreement with the formula obtained by 
Zaremba for the sound velocity of a homogeneous cigar- 
shaped condensate in the TF regime Il28ll . 

In the quasi- ID mean field regime with no vortices 
(Xi < 1, (3 q = 1) Eq. (|50j reduces to 



mc 1D 



hw 



2an\(z) 



1, 



(52) 



which is also the correct result as follows from the substi- 
tution of Eq. (|48j into Eq. (l49l) . 

For the particular case of a homogeneous cigar-shaped 
condensate (no axial confinement and n\ constant) the au- 
thors of Ref. lE3ll have obtained an analytical expression 
for the axial sound velocity that reproduces correctly the 
perturbative result and, to a good approximation (with a 




FIG. 1 : (Color online) Theoretical prediction for the ground-state 
properties (in units of hw z ) of arbitrary condensates with q = 
in different trap geometries A (solid lines). The open circles are 
the exact numerical results. 



relative error less than 3%), also the TF result, interpolat- 
ing between these two regimes. This expression has been 
obtained using an approach different from that used here 
and, in fact, it exhibits a functional dependence on arii that 
is very different from that in Eq. < T50b - However, particular- 
izing Eq. ( TSOb for an axially homogeneous condensate one 
finds that, in this case, both expressions are in quantitative 
agreement within 3.75%. 



IV. NUMERICAL RESULTS 

To verify the predictions of our model we have numeri- 
cally solved the stationary Gross-Pitaevskii equation (TJl by 
using a pseudospectral method evolving in imaginary time. 
Figure Q] shows the ground-state properties of condensates 
with q = and an arbitrary number of particles (xo) in 
different trap geometries (A). The open circles are exact 
numerical results obtained from the Gross-Pitaevskii equa- 
tion. The solid lines are the theoretical predictions (in units 
of hu z ) obtained from Eqs. (fT9ll-(l24l). We have solved 
the polynomial equation (l20l) by using a symbolic software 
package. As already said, this is a trivial computational 
task that only requires one to type in a single instruction. 

As is evident from Fig. [T] the agreement is very good in 
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FIG. 2: (Color online) Theoretical prediction for the ground-state 
properties (in units of hcu z ) of arbitrary condensates with a q — 
1 vortex in different trap geometries A (solid lines). The open 
circles are the exact numerical results. 



all trap geometries (typically better than 1%). In particular, 
although our formalism is cylindrically symmetric it can 
describe accurately the properties of spherical condensates 
(A = 1). For trap anisotropies higher than those considered 
in Fig. [T]one can make use of the analytical solutions (l26l l 
and (l40b found above. These cases will be examined below 

(Fig. H>. 

In Figs. [2H1] we show the ground-state properties of 
condensates containing an axisymmetric vortex of charge 
q = 1, 2, and 4, respectively. As before, the solid lines 
are the theoretical results obtained from Eqs. (fT9t-d24t. 
These figures show that, regardless of the number of par- 
ticles (xo) and trap geometry (A), our model, despite its 
simplicity, can also reproduce very accurately the physical 
properties of condensates containing an axisymmetric vor- 
tex. This is remarkable because, rather crudely, we have 
incorporated the effect of the vortex into a uniform mean 
density n 2 (with a renormalized interaction strength gn^ 1 ) 
localized in the inner cylinder r± < r° . In turn, the radius 
r° follows from Eq. ©, which cannot account for the ef- 
fect of the mean-field interaction energy. However, for con- 
densates of intermediate size, the interaction energy is no 
longer negligible in comparison with the kinetic energy. In 
fact, it plays an important role in determining the size of the 
vortex core, which decreases as N increases. Equation (O 



FIG. 3: (Color online) Same as Fig. |2]for arbitrary condensates 
with a q — 2 vortex. 
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FIG. 4: (Color online) Same as Fig. [2] for arbitrary condensates 
with a q — 4 vortex. 
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FIG. 5: (Color online) (a)-(b): Theoretical prediction for the 
ground-state properties (in units of fko±) of arbitrary cigar- 
shaped condensates with A « 1 and different vortex charges q 
(solid lines). The open circles are exact numerical results ob- 
tained with A = 0.1. (c)-(d): Theoretical prediction for the 
ground-state properties (in units of huj z ) of arbitrary disk-shaped 
condensates with A 3> 2(|g| + 1) and different vortex charges q 
(solid lines). The open symbols are exact numerical results ob- 
tained with A = 100. 



cannot incorporate this correction, which is proportionally 
more important for a high q and an intermediate number 
of particles (for a sufficiently large N the overall contribu- 
tion of the vortex becomes negligible). Thus, one expects 
the formulas above to be less accurate as q increases. This 
can be appreciated from Fig. [4] which shows that for con- 
densates with a q = 4 vortex and a Xo sufficiently large 
that €i nt > ejdn the theoretical predictions are slightly less 
accurate than those corresponding toag=lora(7 = 2 
vortex. 

Next, we consider the ground-state properties of con- 
densates confined in disk-shaped traps satisfying A S> 
2(| q\ + 1) and cigar-shaped traps satisfying A <C 1. As 
already seen, in these limiting cases the theoretical results 
become independent of A and can be directly obtained from 
explicit analytical formulas. In Fig. [5] we show the chem- 
ical potentials and interaction energies obtained from Eqs. 
(|2oTi-(|28Ti and d40b — d42l) . along with exact numerical re- 
sults. The corresponding kinetic and potential energies fol- 
low immediately from Eqs. (l24l . 
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FIG. 6: (Color online) (a): Theoretical prediction for the con- 
densate density per unit length rii (z) of arbitrary xi = 1 cigar- 
shaped condensates with A < 1 and different vortex charges q 
(solid lines). The open symbols are exact numerical results ob- 
tained with A = 0.1. (b): Theoretical prediction for the conden- 
sate density per unit area n 2 (r_[_) of arbitrary \i — 1 disk-shaped 
condensates with A S> 2(|g| + 1) and different vortex charges q 
(solid line). The open symbols are exact numerical results ob- 
tained with A = 100. 



As Fig. |3a) reflects, the chemical potential increases 
with the vortex charge, a consequence of the larger kinetic 
and potential energies associated with multiply quantized 
vortex states. However, as shown in Fig. |3b), the op- 
posite occurs with the mean interaction energy, which de- 
creases as q increases because of the dilution effect that the 
centrifugal barrier produced by the vortex has on the mean 
condensate density. The small errors that can be appreci- 
ated in Fig. [3b) are due, in part, to the fact that, for q ^ 0, 
the approximate solution ( |40t incorporates a residual er- 
ror of order 2-3 %. The exact solution of the polynomial 
equation ( |38l l leads to somewhat better results. 



On the other hand, Figs. |3c) and [5jd) reflect the fact 
that for highly asymmetric trap geometries satisfying A S> 
2{\q\ + 1) the vortex contribution to the condensate prop- 
erties becomes negligible. 

Figure |6ta) shows the condensate density per unit length 
rii(z) of arbitrary cigar-shaped condensates with A <C 1 
[obtained from our analytical formulas (|40|) and d43l)l. 
while Fig. |6jb) shows the condensate density per unit area 
n 2( r ±) °f arbitrary disk-shaped condensates with A S> 
2(\q\ + 1) [obtained from Eqs. <[26) and d29Vl. 

The good agreement between theoretical and exact re- 
sults demonstrates that the formulas derived above are ap- 
plicable for any trap geometry. 
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V. CONCLUSION 

In a previous work we derived very accurate approxi- 
mate analytical expressions for the ground-state properties 
of trapped spherical, cigar-shaped, and disk-shaped con- 
densates with an arbitrary number of atoms in the mean- 
field regime ll26ll . In this work we have extended our pre- 
vious proposal and have derived general approximate for- 
mulas that provide with remarkable accuracy the ground- 
state properties of any mean-field scalar Bose-Einstein 
condensate with short-range repulsive interatomic interac- 
tions, confined in arbitrary cylindrically symmetric har- 
monic traps, and even containing a multiply quantized ax- 
isymmetric vortex. 

In the appropriate limits (corresponding to cigar- shaped 
and disk-shaped condensates) the ground-state properties 
follow from explicit analytical formulas that generalize 
those obtained in Ref. [26]. In the general case, how- 
ever, one has to solve a quintic polynomial equation. In 
this regard, it is important to note that while solving the 
GP equation can be a complex computational problem (es- 
pecially in highly asymmetric trap geometries), solving a 
polynomial equation is a trivial computational task. Us- 
ing the built-in capabilities of symbolic software packages 
such as MATHEMATICA or MATLAB one obtains an instan- 
taneous result after typing in a single instruction. 

The model presented in this work is essentially a con- 
venient approximation method motivated by two simple 
ideas: (i) There exists a direct relation between the number 
of particles and the size of a trapped BEC and (ii) the con- 
tribution from the harmonic oscillator energy to the chem- 
ical potential cannot be smaller than the zero-point energy. 
Applying these simple ideas one finds useful formulas of 
great generality that provide the condensate ground-state 
properties in terms of the correct physically relevant mag- 
nitudes. Moreover, even though no freely adjustable pa- 
rameters are introduced, in all cases the formulas obtained 
reproduce simultaneously with a remarkable accuracy the 
condensate chemical potential and the interaction energy 
and, as a consequence, the kinetic and the potential energy 
as well. And this is true for mean-field condensates with 
any number of particles, confined in any axisymmetric har- 
monic trap, and even containing an axisymmetric vortex. 

Finally, note that the results of this work are applicable 
in general to any nonlinear system characterized by the sta- 
tionary nonlinear Schrodinger equation (TJ). 

This work has been supported by MEC (Spain) and 
FEDER fund (EU) (Contract No. Fis2005-02886). 



T Electronic address: vdelgado@ull.es 
[1] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wie- 

man, and E. A. Cornell, Science 269, 198 (1995). 
[2] K. B. Davis, M. -O. Mewes, M. R. Andrews, N. J. van 

Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Phys. 

Rev. Lett. 75, 3969 (1995). 
[3] C. C. Bradley, C. A. Sackett, and R. G. Hulet, Phys. Rev. 

Lett. 78, 985 (1997). 
[4] For a review see, for example, F. Dalfovo, S. Giorgini, L. P. 

Pitaevskii, and S. Stringari, Rev. Mod. Phys. 71, 463 (1999). 
[5] E. P. Gross, Nuovo Cimento 20, 454 (1961); J. Math. Phys. 

4, 195 (1963); L. P. Pitaevskii, Zh. Eksp. Teor. Fiz. 40, 646 

(1961) [Sov. Phys. JETP 13, 451 (1961)]. 
[6] G. Baym and C. J. Pethick, Phys. Rev. Lett. 76, 6 (1996). 
[7] F. Dalfovo, L. Pitaevskii, and S. Stringari, Phys. Rev. A 54, 

4213 (1996). 

[8] E. Lundh, C. J. Pethick, and H. Smith, Phys. Rev. A 55, 
2126 (1997). 

[9] A. L. Fetter and D. L. Feder, Phys. Rev. A 58, 3185 (1998). 
[10] V. M. Perez-Garcfa, H. Michinel, J. I. Cirac, M. Lewenstein, 

and P. Zoller, Phys. Rev. A 56, 1424 (1997). 
[11] L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A 65, 

043614 (2002); L. Salasnich, Laser Phys. 12, 198 (2002). 
[12] A. L. Fetter, J. Low Temp. Phys. 106, 643 (1997). 
[13] K. K. Das, Phys. Rev. A 66, 053612 (2002). 
[14] P. Schuck and X. Vinas, Phys. Rev. A 61, 43603 (2000). 
[15] E. Timmermans, P. Tommasini, and K. Huang, Phys. Rev. 

A 55, 3645 (1997). 
[16] D. S. Rokhsar, Phys. Rev. Lett. 79, 2164 (1997). 
[17] S. Sinha, Phys. Rev. A 55, 4325 (1997). 
[18] M. Urban, P. Schuck, and X. Vinas, Eur. Phys. J. D 27, 147 

(2003). 

[19] M. Olshanii, Phys. Rev. Lett. 81, 938 (1998). 

[20] D. S. Petrov, G. V. Shlyapnikov, and J. T. M. Walraven, 

Phys. Rev. Lett. 85, 3745 (2000). 
[21] V. Dunjko, V. Lorent, and M. Olshanii, Phys. Rev. Lett. 86, 

5413 (2001). 

[22] C. Menotti and S. Stringari, Phys. Rev. A 66, 043610 
(2002). 

[23] L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A 69, 
045601 (2004). 

[24] D. S. Petrov, M. Holzmann, and G. V. Shlyapnikov, Phys. 

Rev. Lett. 84, 2551 (2000). 
[25] A. Gorlitz, J. M. Vogels, A. E. Leanhardt, C. Raman, T. L. 

Gustavson, J. R. Abo-Shaeer, A. P. Chikkatur, S. Gupta, S. 

Inouye, T. Rosenband, and W. Ketterle, Phys. Rev. Lett. 87, 

130402 (2001). 

[26] A. Munoz Mateo and V. Delgado, Phys. Rev. A 74, 065602 
(2006). 

[27] M. Abramowitz and I. Stegun, Handbook of Mathematical 

Functions (Dover, New York, 1972). 
[28] E. Zaremba, Phys. Rev. A 57, 518 (1998). 



* Electronic address: ammateo@ull.es 



